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Abstract 

Breaking  internal  waves  in  the  vicinity  of  topography  can  reach  heights  of  over  100  m 
and  are  thought  to  enhance  basin-wide  energy  dissipation  and  mixing  in  the  ocean.  The 
scales  at  which  these  waves  are  modelled  often  include  the  breaking  of  large  waves  (10s 
of  meters),  but  not  the  turbulence  dissipation  scales  (centimeters).  Previous  approaches 
to  parameterize  the  turbulence  have  been  to  use  a  universally  large  viscosity,  or  to  use 
mixing  schemes  that  rely  on  Richardson-number  criteria. 

A  simple  alternative  is  presented  that  enhances  mixing  and  viscosity  in  the  pres¬ 
ence  of  breaking  waves  by  assuming  that  dissipation  is  governed  by  the  equivalence  of 
the  density  overturning  scales  to  the  Ozmidov  scale  (e  =  LfN 3,  where  Lj  is  the  size 
of  the  density  overturns,  and  N  the  stratification).  Eddy  diffusivities  and  viscosities  are 
related  to  the  dissipation  by  the  Osborn  relation  ( Kz  =  FeN  2)  to  yield  a  simple  param¬ 
eterization  Kz  =  r LfN,  where  F  ss  0.2  is  the  flux  coefficient.  This  method  is  compared 
to  previous  schemes  for  flow  over  topography  to  show  that,  when  eddy  diffusivity  and 
viscosity  are  assumed  to  be  proportional,  it  dissipates  the  correct  amount  of  energy,  and 
that  the  dissipation  reported  by  the  mixing  scheme  is  consistent  with  energy  losses  in 
the  model.  A  significant  advantage  of  this  scheme  is  that  it  has  no  tunable  parameters, 
apart  from  the  turbulent  Prandtl  number  and  flux  coefficient.  A  disadvantage  is  that  the 
overturning  scales  of  the  turbulence  must  be  relatively  well-resolved. 


1.  Introduction 

There  is  considerable  interest  in  how  topography  interacts  with  stratified  flows  to 
produce  internal  waves  and  turbulence.  In  particular  the  role  of  internal  waves  pro¬ 
duced  by  tides  and  lee  waves  in  flows  over  topography  have  been  examined.  In  general 
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these  have  been  treated  with  large-scale  general  circulation  models  (i.e.  POM  Merri- 
field  and  Holloway  (2002)  or  ROMS),  or  with  specialized  non-hydrostatic  codes  (MIT- 
gcm  Legg  and  Adcroft  (2003),  SUNTANS  Venayagamoorthy  and  Fringer  (2005)).  At 
the  other  end  of  the  spectrum  LES  or  DNS  simulations  have  been  carried  out  on  small 
scales  that  resolve  turbulence. 

Here  we  are  interested  in  the  scales  typical  of  flow  over  deep  ocean  ridges,  like 
Hawaii  (Klymak  et  al.,  2008),  and  at  continental  slopes,  like  Oregon  (Nash  et  ah, 

2007) .  These  flows  are  deep,  up  to  4000  m,  and  exhibit  features  on  a  variety  of 
scales,  from  low-mode  internal  tides  that  span  the  whole  water  column,  to  break¬ 
ing  non-linear  waves  near  abrupt  changes  in  the  topographic  slopes.  These  breaking 
waves  lead  to  density  inversions  that  can  be  over  150-m  tall,  with  dissipation  rates  £  > 
5  x  1 0  s  m2s“3  in  stratifications  N2  ~  1 0  6  s  2.  These  flows  have  turbulent  Reynolds 
numbers  exceeding  106,  and  buoyancy  Reynolds  numbers  Ret,  =  e/vN2  >  104,  and 
Kolmogorov  scales  on  the  order  of  1 0  3  m.  To  capture  the  full  range  of  scales  would 
require  106  grid  cells  in  each  dimension.  Instead,  to  study  these  phenomena,  we  have 
made  extensive  use  of  a  relatively  efficient  class  of  2-D  simulations  that  allow  good 
resolution  in  the  vertical  (O(10m))  and  horizontal  (0(100  m)),  such  that  the  large-scale 
forcing  and  subsequent  breaking  of  internal  waves  can  be  simulated.  Direct  numer¬ 
ical  simulation  methods  are  prohibitively  expensive  for  exploration  of  the  parameter 
spaces  in  which  these  phenomena  are  forced,  therefore  the  turbulence  dissipation  in 
these  simulations  needs  to  be  parameterized. 

Two  approaches  to  parameterizing  turbulence  at  this  scale  of  modelling  have  been 
used.  The  first  is  to  use  a  high  vertical  viscosity  A,  ~  1 0  2  m2s“ 1  (Legg  and  Klymak, 

2008)  or  A,  -  10  1  mV  (Legg  and  Huijts,  2006;  Legg  and  Adcroft,  2003).  This  has 
the  benefit  of  being  simple,  and  yielding  the  turbulence  dissipation  in  the  flow 


As  we  argue  below,  this  scheme  depends  strongly  on  the  choice  of  Az.  Even  if  it  is 
tuned  to  the  breaking  waves,  it  can  exaggerate  dissipation  in  the  midfield  where  the 
shear  is  strong  but  not  strong  enough  to  excite  shear  instability  or  breaking. 

The  alternative  has  been  schemes  that  have  enhanced  mixing  based  on  the  value  of 
the  gradient  Richardson  number  Ri  =  (dlJiJdz.)'  2 A' 2  (where  {//,  is  the  horizontal  com¬ 
ponent  of  velocity,  and  N  the  buoyancy  frequency).  Some  schemes  depend  on  a  crit¬ 
ical  Richardson  number  below  which  the  turbulence  is  increased,  such  as  the  Mellor- 
Yamada  scheme  used  here  (Mellor  and  Yamada,  1982),  while  more  recent  schemes 
remove  the  necessity  for  a  critical  Richardson  number  (Galperin  et  ah,  2007;  Canuto 
et  ah,  2008).  In  all  such  schemes  the  production  rate  of  turbulence  kinetic  energy 
is  assumed  to  be  P  =  Az(dUi,/dz)2  where  A,  is  a  turbulent  vertical  diffusivity  meant 
to  represent  unresolved  eddies.  The  problem  with  these  schemes  in  simulations  with 
resolved  breaking  waves  is  that  the  turbulent  eddies  are  partially  resolved  and  drive 
overturning  so  that  Ri  1  <  0  is  resolved  by  the  model.  All  the  schemes  introduce  an 
arbitrarily  large  viscosity  for  negative  Richardson  number  and  e  calculated  from  the 
local  shear  can  be  unreasonably  high  or  low,  depending  on  this  arbitrary  choice. 

Here  we  present  a  simple  local  scheme  for  mixing  in  breaking  regions  based  on 
the  observed  correlation  between  the  size  of  the  convecting  overturn  and  the  Ozmidov 
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scale.  The  Ozmidov  scale  is  related  to  the  rate  of  turbulence  dissipation  by  e  =  L^N  \ 
Energy  arguments  and  observational  evidence  indicates  that  the  size  of  convectively 
unstable  vertical  displacements  in  a  turbulent  patch  L/.  is  approximately  equivalent  to 
the  Ozmidov  scale:  L,  «  L„  (Dillon,  1982;  Wesson  and  Gregg,  1994;  Mourn,  1996). 
This  brushes  over  significant  changes  in  the  dissipation  during  the  life  of  an  overturn 
(i.e.  Gargett  and  Mourn,  1995;  Smyth  et  al.,  2001),  but  is  a  rough  average  dissipation 
rate.  The  correspondence  between  dissipation  rates  and  the  size  of  overturns  in  convec¬ 
tive  instabilities  like  those  found  over  Hawaii  or  in  fjords  appears  to  be  robust  (Klymak 
and  Gregg,  2004;  Klymak  et  al.,  2008).  We  present  this  scheme  as  a  bridge  between 
large  scale  models  that  do  not  resolve  breaking  waves,  and  small-scale  large-eddy  or 
direct  numerical  simulations. 

Below  we  implement  this  simple  scheme  whereby  the  turbulent  viscosity  Az  is  cal¬ 
culated  from  the  size  of  overturns  driven  by  the  breaking  waves.  We  compare  the  en¬ 
ergy  dissipation  predicted  by  the  parameterized  Az  to  the  energy  lost  from  the  model  for 
the  proposed  scheme  and  compare  to  the  constant-/!,,  runs  and  those  using  the  Mellor- 
Yamada  parameterization,  a  widely  used  Richardson-number  scheme.  Two  idealized 
cases  of  interest  are  considered,  steady  and  oscillating  tidal  flow  over  an  obstacle. 

2.  Methods 

The  model  used  here  is  the  MITgcm  (Marshall  et  al.,  1997;  Legg  and  Klymak, 
2008).  We  use  a  2-dimensional  topography,  with  a  stretched  horizontal  co-ordinate 
system.  For  most  of  the  runs  here  water  depth  H  =  2000  m,  and  vertical  resolution 
was  200  points  with  8 zm  =  10  m;  a  few  runs  were  made  with  H  =  1300  and  1650  m. 
The  horizontal  domain  was  174  km  over  240  grid  cells.  The  inner  80  grid  cells  were 
spaced  100  m  apart,  and  then  the  grid  was  telescoped  linearly  so  that  for  the  outer  cells 
Ax  =  2  km.  The  obstacle  in  all  cases  is  a  Gaussian  shape,  height  from  the  seafloor 
given  by  h  =  hmexp(—x2/W2).  The  width  W  introduces  an  aspect  ratio  to  the  problem 
a o  =  hm/W.  The  model  was  run  in  hydrostatic  mode  for  numerical  efficiency.  Exper¬ 
iments  with  non-hydrostatic  code  did  not  reveal  substantive  differences  in  the  features 
of  interest  here  (see  section  3.2  below  for  a  comparison).  That  is  somewhat  counter¬ 
intuitive,  since  our  proposed  scheme  depends  on  breaking  waves  to  provide  the  tur¬ 
bulence.  However,  the  non-hydrostatic  terms  in  the  vertical  momentum  equation  only 
contribute  to  the  evolution  of  the  breaking,  not  to  its  actual  onset.  This  would  require  a 
considerably  more  isotropic  simulation  grid  than  desirable  for  these  scales,  and  would 
be  amenable  to  a  more  isotropic  mixing  scheme  as  well.  For  most  runs,  horizontal 
viscosities  and  diffusivities  were  kept  as  low  as  numerically  feasible,  at  1 0  4  m2s“1, 
except  where  noted. 

2.1.  Vertical  Turbulence  Schemes 

Constant  Viscosity.  These  runs  compare  with  Legg  and  Huijts  (2006)  and  Legg  and 
Klymak  (2008)  who  used  high  vertical  and  horizontal  viscosities  Az  =  1 0  2  m2s“1, 
and  Ah  =  1 0  1  irr  s  1 .  In  those  papers,  explicit  mixing  was  set  to  zero,  and  handled 
numerically  by  a  Superbee  advection  scheme  (van  Leer,  1979).  In  all  our  simulations 
shown  here,  the  diffusivity  was  set  to  be  the  same  as  the  viscosity,  except  for  a  sen¬ 
sitivity  study  that  shows  small  differences  in  the  dissipation  due  to  the  higher-order 
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scheme  (section  4).  The  constant  viscosity  runs  have  the  simple  advantage  that  dissi¬ 
pation  is  relatively  easy  to  compute  from  the  flow  field  and  local  shears,  and  so  long 
as  the  model  is  well-resolved,  gives  an  accurate  representation  of  the  simulated  dissi¬ 
pation.  The  disadvantage  to  this  scheme  is  that  the  Reynolds  number  can  be  too  low 
to  allow  turbulent  features  to  develop  in  the  first  place,  and  that  it  can  place  too  much 
dissipation  in  regions  that  are  not  expected  to  be  turbulent. 

Mellor  Yamada,  2.0.  The  Mellor-Yamada  formulation  used  by  the  MITgcm  is  a  second- 
order  local  model.  This  version  of  the  scheme  does  not  solve  a  prognostic  equation 
for  the  turbulent  kinetic  energy  (so  this  is  not  to  be  confused  with  what  is  commonly 
called  “Mellor-Yamada  2.5”).  This  scheme  enhances  viscosity  above  an  arbitrary  back¬ 
ground  of  Az  =  10  5irr  s_ 1  if  Ri  =  N2/S 2  <  0.25,  and  viscosity  becomes  very  high  if 
Ri  <  0  in  density  overturns  (figure  1).  With  no  capping,  the  maximum  viscosity  is 
4  x  10  1  m2s_1.  The  implementation  used  here  allows  a  cap  to  this  value,  adding  an 
adjustable  parameter  Amax.  The  production  of  turbulent  kinetic  energy  implied  by  this 
relationship  between  viscosity  and  Richardson  number  is  the  same  as  for  the  more  elab¬ 
orate  MY2.5,  but  in  the  local  scheme  energy  is  assumed  to  be  dissipated  locally  rather 
than  propagating  to  remote  grid  cells;  in  this  paper  we  refer  to  this  as  MY2.0.  While 
we  used  MY2.0  rather  than  MY2.5  because  that  was  the  scheme  already  implemented 
in  the  MITgcm,  it  is  likely  that  similar  results  would  be  found  with  MY2.5  because  the 
production  of  turbulence  is  the  same  in  both  schemes,  and  this  production  term,  rather 
than  the  diffusion  and  advection  of  TKE,  is  the  principal  difference  introduced  by  our 
new  scheme  below.  Similarly  the  choice  of  critical  Richardson  number  has  little  influ¬ 
ence  on  our  comparison,  since  the  flaws  in  the  MY  schemes  which  we  are  addressing 
here  occur  for  Ri  <  0. 


I - 1 - r1 

-10  -5  0 

Ri=N2/S2 


Figure  1:  Mellor-Yamada  2.0  scheme  used  in  the  MITgcm.  Negative  Richardson  num¬ 
bers  mean  that  the  stratification  is  unstable.  The  dashed  line  is  Ri  =  0.25.  The  back¬ 
ground  value  of  A  -  =  10~5  irr  s’  1 . 


New  scheme.  The  high  diffusivity  in  the  MY2.0  scheme  is  appropriate  for  boundary 
layers,  where  a  turbulent  length  scale  can  be  calculated  or  estimated  from  similarity 
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theory.  It  is  not  appropriate  to  turbulence  generated  by  breaking  waves  inside  the  model 
domain.  Instead  we  propose  a  simple  scheme  whereby  each  vertical  profile  in  the 
domain  is  sorted  so  there  are  no  overturns  (Thorpe,  1977)  and  the  distance  each  parcel 
has  been  moved,  A z,  is  related  to  the  dissipation  at  that  location: 

e=(A z)2N3s  (2) 

where  (Vj  is  the  square  of  the  buoyancy  frequency  of  the  sorted  profile,  and  should  be 
greater  than  (or  equal  to)  zero.  The  turbulent  diffusivity  is  then  assumed  to  be 

Kz  =  TeNs2  =  T  (Az)2Ns  (3) 

where  F  =  0.2  is  the  nominal  flux  coefficient  and  K-  has  been  estimated  following 
Osborn  (1980).  The  Osborn  model  is  appropriate  for  turbulent  mixing  in  flows  which 
are  overall  stably  stratified,  as  here.  However,  this  mixing  cannot  take  place  without 
local  overturns,  and  our  focus  here  is  the  scenario  in  which  these  local  overturns  are 
resolved  by  the  model,  but  the  turbulence  is  not.  Because  the  buoyancy  Reynolds 
numbers  are  so  large,  and  in  the  absence  of  compelling  guidance  otherwise,  we  then 
assume  a  turbulent  Prandtl  number  of  one,  so  that  Az=  K-. 

An  example  of  calculating  e  from  an  oscillating  tidal  flow  over  topography  is  shown 
in  figure  2.  Large  isopycnal  displacements  lead  to  non-linear  bores  at  this  location  in 
the  model,  driving  the  turbulence  predicted  by  this  scheme.  An  example  density  profile 
has  overturns  (figure  2b)  that  correspond  to  vertical  displacements  Az  (figure  2c),  yield¬ 
ing  dissipation  rates  of  up  to  e  =  1 0  4  m2  s-3,  in  this  case  at  the  bottom  of  the  overturn 
(figure  2d,  thick  line).  Averaging  these  inhomogeneous  events  over  20  minutes  yields 
an  80-m  thick  region  of  average  dissipations  close  to  1 0  5  m2  s  3  (figure  2d,  thin  line). 
The  flow  regime  and  turbulence  levels  in  this  model  are  similar  to  that  observed  at 
Hawaii  (Klymak  et  ah,  2008). 

There  are  a  number  of  objections  that  could  be  raised  about  this  simple  parame¬ 
terization,  and  the  choices  of  turbulent  Prandtl  number  and  flux  coefficient  could  be 
refined  in  future.  However,  the  goal  here  is  to  move  beyond  the  excessive  diffusivities 
implied  by  Mellor-Yamada-type  schemes  in  overturns,  or  the  arbitrariness  of  constant 
diffusivities.  One  major  limitation  is  that  this  scheme  does  not  work  in  well-mixed 
fluids  where  Ns  =  0,  (we  are  assuming  that  the  flow  is,  following  the  Thorpe  reorder¬ 
ing,  stably  stratified),  a  second  is  that  it  does  not  account  for  shear-instabilities  if  they 
are  not  resolved  in  the  model.  We  have  also  been  somewhat  lax  in  our  definition  of 
the  dissipation  in  an  overturn.  The  proper  definition  of  the  mean  dissipation  over  an 
overturn  would  be  (e)  =  Az2  (N)3,  whereas  our  calculation  will  give  (e)  =  Az2  ((V3). 
It  will  also  distribute  the  dissipation  somewhat  idiosyncratic  ally  where  Az  is  locally 
large.  We  do  not  choose  to  worry  about  these  details  because  the  correlations  that  are 
used  to  justify  the  relation  between  the  Ozmidov  scale  and  the  Thorpe  scale  are  poorly 
constrained  (Mourn,  1996)  and  the  added  “precision”  of  averaging  N  or  e  over  an  over¬ 
turn  requires  significantly  more  computational  complexity  in  identifying  overturning 
patches  and  performing  the  average.  This  is  further  justified  because  the  variation  of 
stratification  in  overturns  tends  to  be  modest  (figure  2b). 

The  scheme  is  implemented  on  the  vertical  profile  from  each  horizontal  grid  cell 
independently.  A  sorted  profile  of  density  is  determined  using  a  simple  insertion-sort 
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Figure  2:  An  example  of  the  proposed  scheme  at  a  single  location  in  a  simulation  of  os¬ 
cillating  stratified  flow  over  an  obstacle.  This  is  taken  atx  =  2.5  km  from  the  simulation 
presented  below  (figure  8).  a)  Dissipation  in  the  water  column,  with  density  contours 
over  top.  Density  contours  are  from  every  50  m  in  the  background  density  profile.  The 
black  vertical  bar  indicates  the  profile  plotted  in  the  other  panels,  b)  density  profile  be¬ 
tween  300  and  450  m,  both  before  and  after  sorting,  c)  vertical  displacement  necessary 
to  make  the  profile  sorted,  d)  Turbulence  dissipation  calculated  for  this  profile  (heavy 
line)  and  the  20-minute  average  around  when  this  profile  was  taken,  e)  Diffusivity  for 
this  profile  (heavy  line)  and  the  20-minute  average. 


algorithm  (Press  et  ah,  1992).  This  algorithm  is  not  well-suited  to  very  turbulent  flows, 
but  for  relatively  laminar  flows  with  rare  overturns  it  should  be  faster  (~  O(n),  where 
n  is  the  number  of  data  points)  than  more  advanced  sorting  algorithms  like  quicksort 
(0(n\ogn))  that  assume  randomized  data.  For  the  simulations  we  are  interested  in 
here,  overturns  are  localized  to  near  the  topography,  so  much  of  the  model  domain 
is  indeed  laminar.  Because  the  sorting  is  local,  the  memory  requirements  are  a  small 
multiple  of  the  number  of  vertical  grid  cells.  For  our  model  runs,  this  sorting  step  had 
an  imperceptible  effect  on  execution  time. 

3.  Examples 

3.1.  Steady  Flow  over  a  Large  Obstacle 

For  the  first  example  we  consider  steady  flow  over  a  large  obstacle,  similar  to  model 
runs  presented  earlier  (Klymak  et  al.,  2010).  Here  a  steady  flow  passes  over  a  fixed 
obstacle  such  that  NH /U0  »  1,  where  H  is  the  water  depth,  and  U0  is  the  steady 
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Figure  3:  Velocities  and  density  contours  for  four  different  turbulence  mixing  schemes 
used  to  simulate  steady  flow  over  large  topography.  All  were  run  with  um/N(H  —  hm)  = 
um/Ndm  =  0.023  a)  Az  =  10"5  m2  s"1,  b)  Az  =  10"2  m2  s"1,  c)  local  Mellor-Yamada, 
and  d)  the  proposed  overturning  scheme. 


velocity.  For  these  runs  the  maximum  viscosity  in  the  Mellor-Yamada  scheme  was 
allowed  to  be  the  maximum  allowed  by  the  default  parameters:  Amax  =  0.4  m2  s"  1 . 

As  discussed  in  Klymak  et  al.  (2010),  this  flow  is  best  parameterized  with  the 
barotropic  flow  over  the  ridge  crest  um  =  UaH / (H  —  hm),  where  hm  is  the  height  of 
the  topography.  For  the  runs  here  N  =  5.2  x  1 0  2  s' ,  //  =  2000  m  ,H  —  hm  =  1000  m, 
and  U0  =  0.03  to  0.18  ms"1.  An  example  of  the  results  for  Ua  =  0.06  ms"1  demon¬ 
strates  the  differences  between  the  mixing  schemes  (figure  3, figure  4).  The  large-scale 
response  is  essentially  the  same  in  all  four  cases,  with  a  stagnant  layer  forming  up-  and 
downstream  of  the  obstacle,  and  acceleration  of  the  flow  above  to  near  um  =  0.06  ms"1 . 
The  high-mode  response  is  relatively  similar  as  well,  with  arrested  waves  just  down¬ 
stream  of  the  crest  that  have  vertical  wavelengths  Xz  sa  27t «  144  m.  There  are, 
however,  differences  in  the  details  of  exactly  how  the  high-mode  response  looks.  The 
low-mixing  run  has  a  number  of  small  scale  instabilities  that  obscure  the  response, 
while  the  other  three  methods  have  differences  in  the  exact  character  of  the  response. 

The  biggest  difference  is  the  dissipation  in  each  run  (figure  4).  We  integrate  the 
dissipation  (figure  4)  for  the  whole  water  column  between  —25  km  <  x  <  25  km,  to  get 
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Figure  4:  Example  dissipation  in  model  runs  with  steady  stratified  flow  over  a  large 
obstacle.  The  totals  are  area-integrated  dissipation  ([m4s~3])  with  DB  the  difference  in 
the  Bernoulli  flux  past  x  =  ±25  km,  and  D  the  reported  dissipation. 


The  dissipation  is  very  small  for  the  Low-mixing  case  (figure  4a)  and  high  for  the 
Mellor-Yamada  case  (figure  4c),  and  similar  for  the  High-mixing  and  New  cases.  The 
greatest  difference  between  the  High  Mix  case  and  the  New  case  is  the  spatial  depen¬ 
dence  of  the  dissipation.  The  dissipation  in  the  High  Mix  case  is  spread  through  the 
domain  whereas  the  dissipation  with  the  New  model  is  localized  in  the  jumps  directly 
in  the  lee  of  the  obstacle  (figure  5).  Given  this,  we  might  expect  that  for  weaker  flows 
the  High-Mix  case  may  start  to  report  higher  dissipation  rates  than  the  New  case. 

The  rate  at  which  energy  is  dissipated  can  be  independently  calculated  from  the 
flow.  These  examples  are  in  steady  state,  so  the  dissipation  of  energy  in  the  flow  is 

DB  =  <j>  Bu  dA  (5) 

where  B  =  \u2  +  P  +  gz  .  The  drop  in  energy  flux  includes  both  explicit  dissipation, 
calculated  as  Z),  and  unknown  numerical  dissipation.  Therefore  we  have  two  types  of 
inaccuracies: 

•  The  model  numerically  dissipates  the  “right”  amount  DB,  but  the  mixing  scheme 
does  not  correctly  report  this  value  ( D  -f-  DB). 

•  The  model  does  not  dissipate  the  “right”  amount  as  predicted  from  the  flow  pa¬ 
rameters  (D b  incorrect). 
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Figure  5:  Dependence  of  the  reported  dissipation  for  the  High  Mix  and  the  New  dissi¬ 
pation  schemes  on  horizontal  location.  The  cumulative  dissipation  shows  how  localized 
the  new  scheme  is. 


Naturally  we  would  like  the  mixing  scheme  to  dissipate  the  “right”  amount  and  “re¬ 
port”  that  amount  accurately.  The  problem  is  that  a-priori  theories  for  the  strength  of 
dissipation  in  any  flow  are  hard  to  come  by,  so  here  we  resort  to  comparing  model 
ensembles.  In  the  examples  above,  we  find  that  the  Low  Mixing  case  (figure  4a)  has  a 
low  model  dissipation  Z)g,  and  that  even  this  low  dissipation  is  grossly  under  predicted 
(D  <  Db).  MY2.0  has  a  slightly  lower  dissipation  than  High  Mix  and  New,  but  it  is 
the  same  order  of  magnitude;  however  it  greatly  over-reports  the  total  dissipation  in  the 
model  (D  >  Du).  Both  High  Mix  and  New  report  similar  dissipations  via  both  methods 
(D  ss  Db). 

A  systematic  comparison  over  a  range  of  resolvable  um/Nhm  demonstrates  the  util¬ 
ity  of  the  various  schemes  (figure  6).  First,  all  the  schemes  dissipate  roughly  the  same 
amount  of  energy  loss  through  the  model  ( Db ,  figure  6a  and  b),  with  the  very  low  mix¬ 
ing  case  ( Az  =  1 0  5  m2s_1)  being  the  prime  exception.  This  loss  of  energy  is  because 
such  low-viscosity  runs  develop  small  scale  numerical  instabilities  that  drain  the  mean 
flow  of  energy.  In  general  the  lower  constant-viscosity  runs  have  slightly  lower  energy 
losses,  but  the  difference  is  not  great.  This  constancy  in  the  energy  lost  is  expected  in 
this  flow  regime  where  the  energy  loss  is  determined  by  the  hydraulic  jump  in  the  lee  of 
the  obstacle,  which  is  itself  determined  by  the  need  of  the  downslope  flow  to  match  the 
downstream  flow.  Unless  the  mixing  scheme  drastically  alters  the  response,  Db  should 
be  a  constant  property  of  the  flow,  and  set  by  matching  the  hydraulic  conditions  at  the 
obstacle  crest  to  the  downstream  conditions  (Baines,  1995). 

The  similarity  of  the  schemes  does  not  persist  when  the  “reported”  dissipation  D  is 
compared  to  the  energy  loss  Db  (figure  6c).  MY2.0  greatly  over-reports  the  dissipation 
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Figure  6:  Comparing  four  different  mixing  schemes,  a)  2-D  integrated  energy  dissipa¬ 
tion  for  different  runs  classified  by  ( Nhm/um)~l ,  as  calculated  by  the  drop  of  energy 
over  the  obstacle  Dg.  b)  The  ratio  of  the  dissipation  Dg  to  D„  =  |jt (H  —  hm)u%t,  which 
is  a  relatively  good  scaling  for  this  flow  regime,  c)  The  ratio  of  the  model  reported 
dissipation  D ,  and  the  observed  energy  drop  Dg.  The  “No-mixing”  and  “MY”  runs 
were  only  made  at  three  flow  speeds. 


(squares).  Most  of  the  dissipation  is  occurring  in  overturns  in  the  lee  of  the  obstacle 
(figure  4c),  so  4  -  «  1 0  1  m2  s  1  which  is  too  high.  This  leads  to  the  energy  lost  from 
the  mean  flow  e  +  K:N2  =A-S 2  to  be  over-reported  because  the  shears  are  significantly 
reduced  over  a  timestep;  the  time-scale  to  diffuse  the  10  m  shears  is  approximately  10 
s,  which  is  less  than  the  12-s  timestep. 

Similarly,  the  constant-viscosity  runs  can  over-  or  under-report  the  dissipation  (fig¬ 
ure  6c,  diamonds).  If  the  viscosity  is  too  low,  numerical  dissipation  makes  up  the 
difference.  If  it  is  too  high,  over-reporting  occurs.  Note  that  Az  =  1 0  2  m2  s'1  does  not 
do  a  bad  job  of  reporting  the  dissipation  (figure  6c,  medium  gray  diamonds),  but  has  a 
distinct  trend  to  under-report  the  more  turbulent  runs  and  over-report  the  less  turbulent 
ones.  For  this  flow  regime  A,  =  1 0  2  m2s_1  is  a  fortuitous  choice  for  the  constant 
viscosity.  The  distribution  of  vertical  viscosities  using  the  New  scheme  indicates  why 
this  is  the  case  (figure  7)  with  the  median  of  the  PDF  very  near  Az  =  1 0  2  m2s_1. 
However,  as  noted  above,  we  do  not  feel  this  puts  the  dissipation  in  the  right  places  in 
the  water  column,  and  it  is  still  an  arbitrary  choice  that  would  not  be  valid  if  the  flow 
regime  were  different.  We  could  improve  the  MY  scheme  with  a  similar  tuning,  but 
again,  it  would  be  tuning  parameters  for  a  particular  regime. 
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Figure  7:  Distribution  functions  of  A-  for  the  New  mixing  scheme  for  the  steady  flow 
over  an  obstacle  case,  with  increasing  flow  velocities. 


3.2.  Oscillating  flow  over  a  Gaussian  bump 

The  original  motivation  for  this  study  was  oscillating  flow  over  a  submarine  ridge. 
Here  we  consider  four  ridges,  hm  =  300,500, 1000  and  1500  —  m  high  in  H  =  2000  m 
of  water.  All  the  ridges  are  Gaussian,  with  W  =  10  km.  Tidal  forcing  every  12.4  h 
on  the  boundaries  induces  U0  =  0.04,0.08,0.12,0.16,0.20,  and  0.24  ms-1  barotropic 
flow  in  the  deep  water,  and  Uj  =  U0H / (H  —  hm)  barotropic  flow  over  the  obstacle. 
All  runs  were  made  with  a  constant  stratification  N  =  5.2  x  1 0  3  s"1,  and  carried  out 
with  the  same  grid  configuration  used  for  the  steady  case.  Here  we  used  a  maximum 
viscosity  for  the  MY  2.0  scheme  of  Amax  =  1 0  2  m2  s'. 

As  above,  energy  is  diagnosed  by  comparing  the  energy  divergence  in  the  model  to 
the  diagnosed  dissipation  (equation  (5)).  In  a  tidal  flow  this  budget  is  never  in  steady 
state,  introducing  another  term: 

DE(t)  =  -jtj  E  dV  +  j Bu  dA,  (6) 

where  E  =  ^u2  +gz  is  the  energy  density.  The  terms  on  the  right  hand  side  are  large, 
but  almost  balance  to  give  a  small  dissipative  residual  D/  ,  so  it  is  imperative  that  these 
terms  are  diagnosed  at  every  time  step  in  the  model,  otherwise  small  errors  in  phase  be¬ 
tween  two  sides  of  the  volume  dominate  actual  divergences  and  good  energy  estimates 
cannot  be  attained.  One  could  linearize  the  energy  sinks  and  sources  (i.e.  by  using  a 
linear  potential  energy  Gill,  1982),  but  this  approach  requires  us  to  assume  the  energy 
content  in  the  volume  is  in  steady  state,  which  is  not  true  due  to  the  background  po¬ 
tential  energy  changing  due  to  mixing.  A  proper  available  potential  energy  calculation 
(Lamb,  2008)  may  yield  better  results,  but  as  that  requires  expensive  global  sorting,  we 
decided  that  it  made  more  sense  to  simply  do  the  absolute  energy  balance  completely. 

It  should  also  be  noted  that  we  do  not  reach  a  complete  energetic  steady  state  with 
these  runs  in  six  tidal  periods  because  slow  high-mode  waves  do  not  reach  the  edge  of 
the  control  volume.  The  runs  could  be  made  longer,  but  at  significant  computational 
cost. 
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Tidal  flow:  U  =0.12  m/s,  h  =  1500  m 
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Figure  8:  Hourly  snapshots  of  energy  dissipation  using  the  New  scheme  for  forced 
semi-diurnal  barotropic  tidal  flow  over  a  Gaussian  bump.  Peak  flow  is  in  the  second 
frame  in  the  positive-x  direction. 


This  issue  with  diagnosing  the  energy  is  demonstrated  in  figure  9,  where  the  dif¬ 
ference  between  a  fully  diagnosed  energy  budget  (a-d)  and  one  from  snapshots  taken 
every  20  minutes  (e-h)  are  compared  for  a  volume  that  is  ±70  km  from  the  obstacle 
crest.  The  energy  fluxes  into  and  out  of  this  boundary  are  very  large  (though  the  ex¬ 
act  magnitude  depends  on  the  definition  of  z  =  0,  and  is  therefore  not  unique).  The 
difference  between  the  energy  divergence  is  balanced  almost  completely  by  the  rate 
of  change  of  energy  in  the  volume.  The  residual  is  very  small  (d)  and  not  properly 
represented  in  the  20-minute  snapshots  (h). 

This  comparison  also  makes  the  point  that  the  high-viscosity  runs  do  not  do  a  very 
good  job  of  tracking  the  energy  dissipation  in  time  (figure  9d).  In  contrast,  the  same 
calculation  made  on  the  new  mixing  scheme  tracks  the  energy  residual  with  time  (fig¬ 
ure  lOd)  more  accurately,  implying  that  the  viscosity  from  the  mixing  scheme  will  yield 
useful  diagnostics  of  energy  loss. 

The  results  over  a  large  number  of  runs  are  similar  to  those  attained  for  the  steady- 
flow  runs  (figure  11).  A  “theoretical”  dissipation,  Dj(H ,hm,U0,N),  is  used  to  nor¬ 
malize  the  dissipations  (and  is  described  in  a  manuscript  in  draft);  while  we  believe 
this  theoretical  dissipation  has  some  usefulness,  the  reader  is  welcome  to  consider  it  a 
scaling  that  happens  to  agree  with  our  New  overturn-based  dissipation.  However,  it  is 
clear  that  using  constant  viscosities  Az  =  1CH2  and  1CT1  m2  s"1  yields  dissipations  that 
depend  on  the  value  of  the  viscosity  (figure  lla-f).  As  in  the  steady  case,  the  reported 
dissipations  ( D ,  figure  lla-c)  have  the  greatest  variation,  whereas  the  energy  loss  in 
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Figure  9:  Time  series  of  energy  budget  for  a  volume  between  x  =  ±70  km  using  con¬ 
stant  Az  -  1 0  2  m2s  *.  a-d)  are  using  a  full  budget  diagnosed  at  every  time  step  of 
the  simulation,  while  e-h)  is  the  same  budget  from  snapshots  taken  every  20  minutes 
(37  times  a  tidal  period),  a)  and  e)  are  the  linear  (dashed)  and  non-linear  (solid)  energy 
fluxes,  the  two  lines  representing  the  in-  and  out-going  fluxes,  b)  and  f)  is  the  total  flux, 
c)  and  g)  are  the  flux  divergences  compared  to  the  rate  of  change  of  energy  in  the  vol¬ 
ume;  these  two  values  lie  almost  one  atop  the  other,  d)  and  h)  is  the  difference  between 
the  energy  flux  divergence  and  the  energy  change  in  the  volume  (solid)  compared  with 
the  prescribed  dissipation  (gray).  Note  the  difference  in  scale  between  d)  and  h). 


the  model  ( De ,  figure  1  ld-f)  tends  to  be  more  constant.  Unlike  the  steady  runs,  the 
MY2.0  scheme  tends  to  under-report  the  dissipation  in  the  oscillating  runs. 
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Figure  10:  As  for  figure  9  except  for  the  overturning-based  dissipation  proposed  here. 
Again,  note  that  d)  and  h)  are  different  scales. 
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Figure  11:  Dissipation  in  tidal  runs  for  different  mixing  schemes:  a-c)  dissipation  re¬ 
ported  by  the  viscocity  D ,  compared  to  a  theoretical  estimate  of  the  dissipation  Dj.  d-f) 
dissipation  calculated  from  an  energy  budget  I)/  ,  g-i)  The  ratio  of  reported  dissipation 
to  diagnosed  dissipation:  D/De- 


The  mixing  scheme  proposed  here  has  the  best  characteristics  of  the  schemes  tested. 
It  underestimates  the  dissipation  slightly  for  the  smallest  (and  thus  least  turbulent)  ridge 
(figure  1  Id),  and  overestimates  slightly  for  the  most  turbulent  (figure  Ilf).  The  size  of 
the  lee  waves  that  drive  much  of  the  turbulence  can  be  estimated  from  the  tidal  flow 
and  stratification  as  An  =  2% Ut/N.  If  there  is  insufficient  resolution  in  the  model  to 
resolve  these  waves  then  there  is  little  hope  they  will  break  and  be  well-represented  by 
the  mixing  scheme. 

Again,  the  differences  in  where  the  dissipation  takes  place  are  the  same  as  for  the 
steady  flow  case:  the  high-viscosity  spreads  the  dissipation  out  over  a  larger  region 
where  tidal  energy  creates  shears  in  the  flow  (figure  12a,b).  These  shears,  however,  are 
stable  to  shear  instabilities  (Ri  >  0.25),  and  therefore  do  not  produce  extra  dissipation 
in  the  MY2.0  (figure  12c,d)  or  the  new  mixing  scheme  (figure  12e,f).  All  three  schemes 
predict  strong  mixing  in  the  lee  of  the  obstacle  during  both  phases  of  the  tide. 

As  a  final  note,  we  tested  the  importance  of  hydrostaticity  for  the  New  mixing 
scheme  in  these  simulations  (figure  13).  For  both  reported  dissipation  and  the  energy 
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Figure  12:  Mean  dissipation  over  a  tidal  period  for  Ut  =  0.16  m  s'1,  h,„  =  1500  m,  and 
N  =  5.2  x  1 0  1  s  ' 1  for  three  mixing  schemes.  Upper  panels  (a,c,e)  are  depth  integrals 
of  the  dissipation  as  a  function  of  distance  from  the  obstacle. 


divergence,  the  response  was  almost  identical  for  the  model  run  in  hydrostatic  and 
non-hydrostatic  configurations.  The  differences  did  seem  to  be  larger  for  larger  forcing 
(figure  13e,f),  perhaps  indicating  that  the  dynamics  inside  the  larger  breaking  waves 
were  playing  a  role,  however,  the  time-smoothed  averages  were  not  significantly  dif¬ 
ferent.  Note  again,  that  the  model  has  a  10:1  aspect  ratio  in  the  grid  cells  near  the 
obstacle,  so  the  lack  of  strong  non-hydrostatic  effects  is  not  surprising.  Experiments 
with  more  horizontal  resolution  may  turn  up  more  substantive  differences,  but  again, 
such  simulations  could  be  made  with  a  more  isotropic  mixing  model  that  acts  on  the 
turbulent  scales  themselves. 

4.  Resolution  and  Advection-Scheme  Limitations 

The  MITgcm,  or  any  numerical  model  that  we  know  of,  does  not  explicitly  try  to 
conserve  mechanical  energy,  so  the  ability  to  attain  near-balances  as  we  have  attempted 
to  show  above,  is  encouraging.  There  are  a  number  of  limitations  to  the  new  scheme 
(and  the  other  schemes)  that  should  be  accounted  for,  and  we  attempt  to  discuss  them 
here. 

First,  the  energy  dissipation  depends  somewhat  on  the  advection  scheme  used. 
Most  of  the  runs  in  this  paper  were  made  with  a  simple  and  efficient  second-order 
scheme.  This  scheme  is  quite  noisy,  and  some  of  this  noise  goes  into  producing  extra 
dissipation  in  the  model.  This  is  seen  most  clearly  when  considering  the  resolution  de¬ 
pendence  of  all  the  schemes  (figure  14).  Here  the  same  runs  were  made  with  finer  res¬ 
olutions.  All  the  schemes  demonstrate  more  dissipation  as  the  resolution  is  increased. 
However,  the  New  scheme  does  so  more  precipitously.  A  2.5-m  resolution  run  shows 
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Figure  13:  Comparison  of  hydrostatic  and  non-hydrostatic  dissipation  and  energy  loss 
from  three  forcings  over  a  Gaussian  obstacle. 


even  more  dissipation,  indicating  that  increasing  vertical  resolution  makes  the  dissipa¬ 
tion  grow  too  much.  It  is  relatively  difficult  to  tell  what  causes  the  extra  dissipation, 
but  it  is  not  just  a  simple  reporting  error  as  it  shows  up  in  both  D  and  /)/,.  We  sug¬ 
gest  that  the  extra  dissipation  is  occurring  because  finer  resolution  allows  more  noise 
to  develop  due  to  inaccuracies  in  the  advection  scheme,  which  makes  more  overturns 
and  therefore  more  dissipation.  The  constant-viscosity  simulations  also  have  increased 
dissipation  with  resolution,  but  the  increase  is  smaller. 

If  we  use  a  nonlinear  limited  advection  scheme  (figure  14b),  in  this  case  the  Su¬ 
perbee  flux-limited  scheme,  we  find  that  the  increase  in  dissipation  with  increased  res¬ 
olution  is  almost  the  same  as  the  constant-viscosity  runs.  However,  the  solutions  are 
smoother  and  the  dissipation  driven  by  overturns  is  now  somewhat  lower  than  before 
for  the  base  10-m  runs. 

We  consider  the  difference  in  dissipations  in  the  New  scheme  more  systematically 
with  a  set  of  experiments  of  tidal  flow  over  a  ridge  using  the  nominal  10-m  resolution 
(figure  15).  In  general,  the  second-order  scheme  has  higher  dissipations  than  the  flux- 
limiting  scheme  for  all  dissipations,  with  the  second  order  scheme  over-reporting  the 
dissipation  (figure  15a,  circles)  and  the  Superbee  scheme  under-reporting  (figure  15a, 
diamonds).  The  two  schemes  tend  to  dissipate  the  same  amount  of  energy  from  the 
model  (De,  figure  15b)  except  for  the  least-turbulent  flows,  where  the  Superbee  scheme 
removes  more  energy  from  the  flow  than  the  2-nd  order.  This  is  consistent  with  the 
Superbee  scheme  having  some  additional  numerical  diffusion. 

As  mentioned  above,  the  mixing  scheme  and  the  simulations  have  a  natural  limit 
where  the  vertical  resolution  becomes  too  coarse  to  resolve  the  waves  in  the  lee  of  the 
obstacle. 
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Figure  14:  Comparison  of  dissipations  at  different  vertical  resolutions  with  a)  2-nd 
order  centered  advection  scheme,  and  b)  Superbee  advection  scheme.  The  triangle  in 
each  panel  is  for  2.5-m  vertical  resolution. 


5.  Summary 

We  have  presented  a  simple  mixing  scheme  for  models  that  resolve  breaking  inter¬ 
nal  waves  based  on  the  Ozmidov  scale  implied  by  these  breaking  waves.  We  expect  this 
scheme  will  be  the  most  useful  for  parameterizing  mixing  in  strong  flows  over  topog¬ 
raphy,  with  particular  application  to  internal  tide  generation  and  reflection  problems.  It 
may  also  be  useful  for  mean-flows  over  rough  topography  such  as  might  occur  in  the 
Antarctic  circumpolar  current. 

There  is  a  slight  cost  to  the  scheme  in  terms  of  memory  and  processor  time.  The 
sorting  algorithm  should  be  relatively  cheap  except  for  very  turbulent  flows. 

We  have  used  some  simplifying  assumptions  that  are  common  in  the  literature  such 
as  that  the  Prandtl  number  is  one  and  that  the  flux  coefficient  F  is  a  constant.  There 
is  significant  interest  in  if  these  assumptions  are  valid  (i.e.  Gargett  and  Mourn,  1995; 
Smyth  et  ah,  2005,  2001)  that  is  beyond  the  scope  of  our  study  here. 

The  major  limitation  of  our  new  mixing  scheme  is  that  the  largest  mixing  events 
should  be  driven  by  breaking  waves  rather  than  unresolved  shear  instabilities.  For 
instance,  in  the  tidal  mixing  over  the  Knight  Inlet  sill  there  is  clear  evidence  of  both 
shear  instability  and  breaking  in  the  lee  wave  (Farmer  and  Armi,  1999;  Klymak  and 
Gregg,  2004).  For  the  flows  discussed  in  this  paper  resolved  shear  instability  was  very 
rare;  the  Mellor-Yamada  mixing  was  almost  always  maximized,  indicating  that  it  was 
responding  to  convective  instability,  rather  than  high  shears  (figure  12).  However,  if 
a  flow  was  being  considered  where  shear  instability  was  also  thought  to  be  important, 
such  as  larger-scale  overflows  (Legg  et  al.,  2009),  it  would  be  very  easy  to  make  a 
hybrid  mixing  scheme  that  also  enhanced  mixing  in  shear  layers.  Rather  than  the  profile 
in  figure  1,  one  could  imagine  the  same  profile  for  Ri  >  0  and  then  an  abrupt  transition 
to  the  overturning  mixing  for  negative  Richardson  Numbers.  Alternate  shear-driven 
schemes,  such  as  that  proposed  by  Jackson  et  al.  (2008),  could  also  be  added. 

The  authors  thank  Alistair  Adcroft  for  helpful  comments  on  a  draft  of  this  pa¬ 
per,  and  two  anonymous  reviewers  and  Boris  Galperin  for  their  very  helpful  com¬ 
ments  on  the  submitted  version.  JK  was  supported  by  ONR  grants  N000 14-08-1- 
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Figure  15:  Dependence  of  dissipation  on  advection  scheme  for  10-m  resolution  os¬ 
cillating  flow  over  an  obstacle:  a)  dissipation  compared  with  energy  divergence  for 
2nd-order  scheme  (circles)  and  superbee  advection  scheme  (diamonds).  The  superbee 
scheme  tends  to  have  some  of  its  dissipation  go  unreported,  the  2-nd  order  scheme  over¬ 
reports  its  dissipation,  b)  the  ratio  of  energy  divergences  in  the  models.  Both  schemes 
lose  approximately  the  same  amount  of  energy  except  at  low  dissipations  where  the 
superbee  scheme  dissipates  more,  c)  Because  it  under-reports,  the  superbee  scheme 
has  less  dissipation  than  the  2-nd  order  scheme. 


0376  and  N00014-08-1-0274  and  NSERC  grant  327920.  SL  was  supported  by  award 
NA08OAR4320752  from  the  National  Oceanic  and  Atmospheric  Administration,  U.S. 
Department  of  Commerce.  The  statements,  findings,  conclusions,  and  recommenda¬ 
tions  are  those  of  the  author(s)  and  do  not  necessarily  reflect  the  views  of  the  National 
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